Fit main model (see exploratory scripts and model comparison), visualize results.
library(tidyverse); theme_set(theme_light(base_size = 12))
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(viridis)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(gganimate)
library(gifski)
library(latex2exp)
library(patchwork)
library(png)
library(RCurl)
library(wesanderson)
library(qwraps2)
library(ggcorrplot)
library(ggridges)
# remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)
# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = "R/analysis/condition_model_cf_cache/html")
# Specify map ranges
ymin = 52; ymax = 60; xmin = 10; xmax = 24
map_data <- rnaturalearth::ne_countries(
scale = "medium",
returnclass = "sf", continent = "europe")
# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
st_crop(map_data,
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))
# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)
# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 8),
legend.text = element_text(size = 8),
legend.title = element_text(size = 8),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin = margin(-5, -5, -5, -5),
strip.text = element_text(size = 8, colour = 'gray10', margin = margin()),
strip.background = element_rect(fill = "grey95")
)
}
# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 6),
strip.text = element_text(size = 8, colour = 'gray10', margin = margin()),
strip.background = element_rect(fill = "gray95"),
legend.position = c(0.7, 0.02),
legend.direction = "horizontal"
)
}
# Make default base map plot
xmin2 <- 303000; xmax2 <- 916000; xrange <- xmax2 - xmin2
ymin2 <- 5980000; ymax2 <- 6450000; yrange <- ymax2 - ymin2
plot_map_raster <-
ggplot(swe_coast_proj) +
xlim(xmin2, xmax2) +
ylim(ymin2, ymax2) +
labs(x = "Longitude", y = "Latitude") +
geom_sf(size = 0.3) +
theme_plot() +
NULL
plot_map_raster_labels <-
plot_map_raster +
annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = 1.9) +
annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = 1.9, angle = 75) +
annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = 1.9) +
annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = 1.9) +
annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = 1.9) +
annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = 1.9, angle = 75) +
annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = 1.9, angle = 75)
# Read the split files and join them
d1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cond_(1_2).csv")
d2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cond_(2_2).csv")
d <- bind_rows(d1, d2)
unique(is.na(d$density_cod))
#> [1] FALSE
unique(is.na(d$density_cod_rec))
#> [1] FALSE TRUE
# Calculate standardized variables
d <- d %>%
mutate(ln_length_cm = log(length_cm),
ln_weight_g = log(weight_g),
year_ct = year - mean(year),
biomass_her_sc = biomass_her,
biomass_her_sd_sc = biomass_her_sd,
biomass_spr_sc = biomass_spr,
biomass_spr_sd_sc = biomass_spr_sd,
density_cod_sc = density_cod,
density_cod_rec_sc = density_cod_rec,
density_fle_sc = density_fle,
density_fle_rec_sc = density_fle_rec,
density_saduria_sc = density_saduria,
density_saduria_rec_sc = density_saduria_rec,
depth_sc = depth,
depth_rec_sc = depth_rec,
depth_sd_sc = depth_sd,
oxy_sc = oxy,
oxy_rec_sc = oxy_rec,
oxy_sd_sc = oxy_sd,
temp_sc = temp,
temp_rec_sc = temp_rec,
temp_sd_sc = temp_sd,
year_f = as.factor(year)) %>%
mutate_at(c("biomass_her_sc", "biomass_her_sd_sc", "biomass_spr_sc", "biomass_spr_sd_sc",
"density_cod_sc", "density_cod_rec_sc",
"density_fle_sc", "density_fle_rec_sc",
"density_saduria_sc", "density_saduria_rec_sc",
"depth_sc", "depth_rec_sc", "depth_sd_sc",
"oxy_sc", "oxy_rec_sc", "oxy_sd_sc",
"temp_sc", "temp_rec_sc", "temp_sd_sc"
),
~(scale(.) %>% as.vector)) %>%
mutate(year = as.integer(year))
unique(is.na(d))
#> weight_g length_cm year lat lon ices_rect sub_div depth oxy temp
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [5,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [6,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [7,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> X Y density_cod density_fle density_cod_data density_fle_data
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE TRUE TRUE
#> [5,] FALSE FALSE FALSE FALSE TRUE TRUE
#> [6,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [7,] FALSE FALSE FALSE FALSE TRUE TRUE
#> density_saduria density_saduria_rec density_cod_rec density_fle_rec
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE
#> [5,] FALSE TRUE TRUE TRUE
#> [6,] FALSE TRUE TRUE TRUE
#> [7,] FALSE FALSE FALSE FALSE
#> depth_rec oxy_rec temp_rec biomass_spr biomass_her density_saduria_sd
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE TRUE TRUE FALSE
#> [3,] FALSE FALSE FALSE TRUE TRUE FALSE
#> [4,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [5,] TRUE TRUE TRUE TRUE TRUE FALSE
#> [6,] TRUE TRUE TRUE TRUE TRUE FALSE
#> [7,] FALSE FALSE FALSE TRUE TRUE FALSE
#> density_cod_sd density_fle_sd depth_sd oxy_sd temp_sd biomass_spr_sd
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE FALSE TRUE
#> [4,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [5,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [6,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [7,] FALSE FALSE FALSE FALSE FALSE FALSE
#> biomass_her_sd ln_length_cm ln_weight_g year_ct biomass_her_sc
#> [1,] FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE TRUE
#> [3,] TRUE FALSE FALSE FALSE TRUE
#> [4,] FALSE FALSE FALSE FALSE FALSE
#> [5,] FALSE FALSE FALSE FALSE TRUE
#> [6,] FALSE FALSE FALSE FALSE TRUE
#> [7,] FALSE FALSE FALSE FALSE TRUE
#> biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc density_cod_sc
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE TRUE FALSE FALSE
#> [3,] TRUE TRUE TRUE FALSE
#> [4,] FALSE FALSE FALSE FALSE
#> [5,] FALSE TRUE FALSE FALSE
#> [6,] FALSE TRUE FALSE FALSE
#> [7,] FALSE TRUE FALSE FALSE
#> density_cod_rec_sc density_fle_sc density_fle_rec_sc density_saduria_sc
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE
#> [5,] TRUE FALSE TRUE FALSE
#> [6,] TRUE FALSE TRUE FALSE
#> [7,] FALSE FALSE FALSE FALSE
#> density_saduria_rec_sc depth_sc depth_rec_sc depth_sd_sc oxy_sc oxy_rec_sc
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [5,] TRUE FALSE TRUE FALSE FALSE TRUE
#> [6,] TRUE FALSE TRUE FALSE FALSE TRUE
#> [7,] FALSE FALSE FALSE FALSE FALSE FALSE
#> oxy_sd_sc temp_sc temp_rec_sc temp_sd_sc year_f
#> [1,] FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE FALSE
#> [5,] FALSE FALSE TRUE FALSE FALSE
#> [6,] FALSE FALSE TRUE FALSE FALSE
#> [7,] FALSE FALSE FALSE FALSE FALSE
unique(is.na(d$density_cod_rec))
#> [1] FALSE TRUE
unique(is.na(d$density_cod))
#> [1] FALSE
# Drop NA covariates for the correlation plot
d <- d %>% drop_na(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc,
density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc)
# Plot correlation between variables
d_cor <- d %>% dplyr::select("biomass_her_sc", "biomass_her_sd_sc",
"biomass_spr_sc", "biomass_spr_sd_sc",
"density_cod_sc", "density_cod_rec_sc",
"density_fle_sc", "density_fle_rec_sc",
"density_saduria_sc", "density_saduria_rec_sc",
"depth_sc", "depth_rec_sc",
"oxy_sc", "oxy_rec_sc",
"temp_sc", "temp_rec_sc")
corr <- round(cor(d_cor), 1)
ggcorrplot(corr, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.3))
ggsave("figures/supp/condition/correlation_variables.png", width = 6.5, height = 6.5, dpi = 600)
pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")
pred_grid <- bind_rows(pred_grid1, pred_grid2)
unique(is.na(pred_grid$density_cod))
#> [1] FALSE
unique(is.na(pred_grid$density_cod_rec))
#> [1] FALSE
pred_grid <- pred_grid %>% mutate(year = as.integer(year),
year_f = as.factor(year))
# Scale the variables with respect to data! First calculate means in data
data_means <- d %>%
dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
density_cod, density_cod_rec,
density_fle, density_fle_rec,
density_saduria, density_saduria_rec,
depth, depth_rec, depth_sd,
oxy, oxy_rec, oxy_sd,
temp, temp_rec, temp_sd) %>%
mutate_all(~mean(.)) %>%
distinct(.keep_all = TRUE)
data_stdev <- d %>%
dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
density_cod, density_cod_rec,
density_fle, density_fle_rec,
density_saduria, density_saduria_rec,
depth, depth_rec,
oxy, oxy_rec,
temp, temp_rec
) %>%
mutate_all(~sd(.)) %>%
distinct(.keep_all = TRUE)
# Before actually scaling, replace the ices-rectangle pelagic values that are NA with the mean across rectangles
# in the sub division. Replace the sub-division pelagic values that are NA with the mean in the year.
# Note that the sub-division values are not the sum of the rectangles (due to missing rectangles), so
# I need to calculate a sub-division mean across rectangles within each sub-division
pred_grid <- pred_grid %>%
group_by(year) %>%
mutate(median_biomass_sprat_across_sub_div = median(biomass_spr_sd, na.rm = TRUE),
median_biomass_herring_across_sub_div = median(biomass_her_sd, na.rm = TRUE)) %>%
ungroup() %>% # Replace sub_divsion NA's with the median across sub_divisions in that year
mutate(biomass_spr_sd = ifelse(is.na(biomass_spr_sd) == TRUE, median_biomass_sprat_across_sub_div, biomass_spr_sd),
biomass_her_sd = ifelse(is.na(biomass_her_sd) == TRUE, median_biomass_herring_across_sub_div, biomass_her_sd)) %>%
group_by(year, sub_div) %>%
mutate(median_biomass_sprat_across_rect_in_sub_div = median(biomass_spr, na.rm = TRUE),
median_biomass_herring_across_rect_in_sub_div = median(biomass_her, na.rm = TRUE)) %>%
ungroup() %>% # Replace rectangle NA's with the median across rectangles in that year and sub-division
mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, median_biomass_sprat_across_rect_in_sub_div, biomass_spr),
biomass_her = ifelse(is.na(biomass_her) == TRUE, median_biomass_herring_across_rect_in_sub_div, biomass_her)) %>%
# Since I still have some NAs (some sub-divisions do not have a single rectangle in some years), I will will those rectangles
# with the sub-division value divided by the number of rectangles in that sub division
group_by(year, sub_div) %>%
mutate(n_rect = length(unique(ices_rect))) %>%
ungroup() %>%
mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, biomass_spr_sd/n_rect, biomass_spr),
biomass_her = ifelse(is.na(biomass_her) == TRUE, biomass_her_sd/n_rect, biomass_her))
pred_grid <- pred_grid %>%
mutate(ln_length_cm = log(1),
year_ct = year - mean(year),
biomass_her_sc = (biomass_her - data_means$biomass_her) / data_stdev$biomass_her,
biomass_her_sd_sc = (biomass_her_sd - data_means$biomass_her_sd) / data_stdev$biomass_her_sd,
biomass_spr_sc = (biomass_spr - data_means$biomass_spr) / data_stdev$biomass_spr,
biomass_spr_sd_sc = (biomass_spr_sd - data_means$biomass_spr_sd) / data_stdev$biomass_spr_sd,
density_cod_sc = (density_cod - data_means$density_cod) / data_stdev$density_cod,
density_cod_rec_sc = (density_cod_rec - data_means$density_cod_rec) / data_stdev$density_cod_rec,
density_fle_sc = (density_fle - data_means$density_fle) / data_stdev$density_fle,
density_fle_rec_sc = (density_fle_rec - data_means$density_fle_rec) / data_stdev$density_fle_rec,
density_saduria_sc = (density_saduria - data_means$density_saduria) / data_stdev$density_saduria,
density_saduria_rec_sc = (density_saduria_rec - data_means$density_saduria_rec) / data_stdev$density_saduria_rec,
depth_sc = (depth - data_means$depth) / data_stdev$depth,
depth_rec_sc = (depth_rec - data_means$depth_rec) / data_stdev$depth_rec,
oxy_sc = (oxy - data_means$oxy) / data_stdev$oxy,
oxy_rec_sc = (oxy_rec - data_means$oxy_rec) / data_stdev$oxy_rec,
temp_sc = (temp - data_means$temp) / data_stdev$temp,
temp_rec_sc = (temp_rec - data_means$temp_rec) / data_stdev$temp_rec,
)
# Plot on map:
ggplot(pred_grid2, aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)
pred_grid %>%
drop_na() %>%
ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)
pred_grid %>%
drop_na(biomass_spr_sd) %>%
ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)
unique(is.na(pred_grid))
#> X Y depth year oxy temp lon lat sub_div ices_rect
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> density_saduria biomass_spr biomass_her biomass_spr_sd biomass_her_sd
#> [1,] FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE
#> density_cod density_fle depth_rec temp_rec oxy_rec density_cod_rec
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> density_fle_rec density_saduria_rec depth_sd temp_sd oxy_sd density_cod_sd
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> density_fle_sd density_saduria_sd year_f
#> [1,] FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE
#> median_biomass_sprat_across_sub_div median_biomass_herring_across_sub_div
#> [1,] FALSE FALSE
#> [2,] FALSE FALSE
#> median_biomass_sprat_across_rect_in_sub_div
#> [1,] FALSE
#> [2,] TRUE
#> median_biomass_herring_across_rect_in_sub_div n_rect ln_length_cm year_ct
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] TRUE FALSE FALSE FALSE
#> biomass_her_sc biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> density_cod_sc density_cod_rec_sc density_fle_sc density_fle_rec_sc
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> density_saduria_sc density_saduria_rec_sc depth_sc depth_rec_sc oxy_sc
#> [1,] FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE
#> oxy_rec_sc temp_sc temp_rec_sc
#> [1,] FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE
pred_grid %>%
drop_na(biomass_spr_sc, biomass_spr_sd_sc, biomass_her_sc, biomass_her_sd_sc) %>%
ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)
spde <- make_mesh(d, xy_cols = c("X", "Y"),
n_knots = 100,
type = "kmeans", seed = 42)
# Plot and save spde
png(file = "figures/supp/condition/spde.png", units = "in", width = 6.5, height = 6.5, res = 300)
plot(spde)
dev.off()
ptm <- proc.time()
mnull <- sdmTMB(
formula = ln_weight_g ~ ln_length_cm,
data = d,
time = NULL,
mesh = spde,
family = student(link = "identity", df = 5),
spatiotemporal = "off",
spatial = "off",
silent = TRUE,
reml = FALSE
)
proc.time() - ptm
#> user system elapsed
#> 1.945 0.263 2.765
# Extract average a and b across the total dataset
a <- tidy(mnull)[1, "estimate"]
b <- tidy(mnull)[2, "estimate"]
a
#> [1] -4.638063
b
#> [1] 2.983866
d$weight_g_avg_pred <- exp(a)*d$length_cm^b
plot(d$weight_g ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$weight_g); abline(a = 0, b = 1, col = "red")
d$le_cren <- d$weight_g / d$weight_g_avg_pred
d$log_le_cren <- log(d$le_cren)
hist(d$le_cren)
mnull_oxy_bp <- sdmTMB(
formula = log_le_cren ~ -1 + year_f + breakpt(oxy_sc) + oxy_rec_sc,
time_varying = NULL,
data = d,
time = "year",
mesh = spde,
family = student(link = "identity", df = 5),
spatiotemporal = "iid",
spatial = "on",
silent = TRUE,
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
mnull_oxy <- sdmTMB(
formula = log_le_cren ~ -1 + year_f + oxy_sc + oxy_rec_sc,
time_varying = NULL,
data = d,
time = "year",
mesh = spde,
family = student(link = "identity", df = 5),
spatiotemporal = "iid",
spatial = "on",
silent = TRUE,
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
AIC(mnull_oxy_bp, mnull_oxy)
min(AIC(mnull_oxy_bp), AIC(mnull_oxy))
#> [1] -144518.8
sanity(mnull_oxy_bp)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
sanity(mnull_oxy)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mnull_oxy)
tidy(mnull_oxy_bp)
nd_bp_simple <- data.frame(
oxy_sc = seq(min(d$oxy_sc), max(d$oxy_sc), length.out = 100),
year = 2019L,
year_f = as.factor(2019),
oxy_rec_sc = 0
)
p_bp_simple <- predict(mnull_oxy_bp, newdata = nd_bp_simple,
se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y"))
ggplot(p_bp_simple, aes(oxy_sc, est,
ymin = est - 1.96 * est_se,
ymax = est + 1.96 * est_se)) +
geom_line() + geom_ribbon(alpha = 0.4)
ggplot(d, aes(oxy_sc, oxy)) +
geom_point() +
geom_hline(yintercept = (1.0700*sd(d$oxy)) + mean(d$oxy)) +
geom_vline(xintercept = 1.0700)
ptm <- proc.time()
mfull <- sdmTMB(
formula = log_le_cren ~ -1 + year_f + biomass_her_sc + biomass_her_sd_sc +
biomass_spr_sc + biomass_spr_sd_sc +
density_cod_sc + density_cod_rec_sc +
density_fle_sc + density_fle_rec_sc +
density_saduria_sc + density_saduria_rec_sc +
depth_sc + depth_rec_sc +
oxy_sc + oxy_rec_sc +
temp_sc + temp_rec_sc,
data = d,
time = "year",
mesh = spde,
family = student(link = "identity", df = 5),
spatiotemporal = "iid",
spatial = "on",
silent = TRUE,
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1)
)
proc.time() - ptm
#> user system elapsed
#> 721.928 20.345 747.627
tidy(mfull, conf.int = TRUE) %>% arrange(desc(abs(estimate)))
tidy(mfull, effects = "ran_par", conf.int = TRUE)
# Check model convergence
sanity(mfull)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
mcmc_res <- residuals(mfull, type = "mle-mcmc", mcmc_iter = 5000, mcmc_warmup = 2000)
png(file = "figures/supp/condition/qq.png", units = "in", width = 6.5, height = 6.5, res = 300)
qqnorm(mcmc_res);qqline(mcmc_res)
dev.off()
d$residuals <- mcmc_res[, 1]
b_fe <- tidy(mfull)
b_fe <- stats::setNames(b_fe$estimate, b_fe$term)
X <- mfull$tmb_data$X_ij
X <- matrix(unlist(X), ncol = length(b_fe))
VarF <- var(as.vector(b_fe %*% t(X))) # variance from fixed-effects
b_ran <- tidy(mfull, "ran_par")
sigma_O <- b_ran$estimate[b_ran$term == "sigma_O"] # spatial standard deviation
sigma_E <- b_ran$estimate[b_ran$term == "sigma_E"] # spatiotemporal standard deviation
VarO <- sigma_O^2 # spatial variance
VarE <- sigma_E^2 # spatiotemporal variance
# Calculate obs - pred
d_r2 <- d
d_r2$pred <- predict(mfull)
d_r2$resid_manual <- d_r2$le_cren - d_r2$pred$est
#hist(d_r2$resid_manual)
VarR <- var(as.vector(d_r2$resid_manual)) # "residual" variance
# Marginal R2s
denominator <- VarF + VarO + VarE + VarR
marg <- VarF/denominator
out <- data.frame(
marginal = marg,
partial_spatial = VarO/denominator,
partial_spatiotemporal = VarE/denominator,
conditional = (VarF + VarO + VarE)/denominator,
all_random = (VarO + VarE)/denominator,
marginal_random_ratio = marg / ((VarO + VarE)/denominator),
random_marginal_ratio = ((VarO + VarE)/denominator) / marg
)
out
# What if I skip the year effect?
b_fe2 <- tidy(mfull) %>% filter(!grepl('year', term))
#> filter: removed 27 rows (63%), 16 rows remaining
b_fe2 <- stats::setNames(b_fe2$estimate, b_fe2$term)
X2 <- mfull$tmb_data$X_ij
X2 <- matrix(unlist(X2), ncol = length(b_fe))[, 28:43] # skip year columns
VarF2 <- var(as.vector(b_fe2 %*% t(X2))) # variance from fixed-effects
denominator3 <- VarF2 + VarO + VarE + VarR
marg3 <- VarF2/denominator3
marg3
#> [1] 0.05794511
# Residuals on map
plot_map_raster +
geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = residuals)) +
scale_colour_gradient2() +
facet_wrap(~ year, ncol = 5) +
labs(color = "residuals") +
NULL
ggsave("figures/supp/condition/spatial_resid.png", width = 6.5, height = 6.5, dpi = 600)
# Residuals vs length and year
ggplot(d, aes(x = length_cm, y = residuals, color = residuals)) +
geom_point(size = 0.1, alpha = 1) +
geom_hline(yintercept = 0, size = 0.5, color = "black", linetype = 2, alpha = 0.5) +
labs(x = "length [cm]", y = "Residuals", color = "") +
stat_smooth(color = "black", size = 0.5) +
facet_wrap(~year, ncol = 5) +
scale_color_gradient2() +
theme_facet_map()
ggsave("figures/supp/condition/residuals_vs_length_and_year.png", width = 6.5, height = 6.5, dpi = 600)
ggplot(d, aes(year, length_cm, z = residuals)) +
stat_summary_2d(bins = 40) +
scale_fill_gradient2()
# ggsave("figures/supp/condition/residuals_vs_length_and_year2.png", width = 6.5, height = 6.5, dpi = 600)
get_index_sims method to avoid slow
bias correctionncells <- filter(pred_grid, year == max(pred_grid$year)) %>% nrow()
#> filter: removed 240,214 rows (96%), 9,239 rows remaining
pred_avg_sim <- predict(mfull, newdata = pred_grid, nsim = 100)
index_avg_sim <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_sims <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)),
return_sims = TRUE) # This is for calculation on samples!
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
write.csv(index_avg_sim, "output/index_avg_sim.csv")
Evaluate sensitivity to using a pred grid with a subset excluding low density areas
# Now do the same but excluding the areas not sampled in the pred grid
ncells_sub <- filter(pred_grid, year == max(pred_grid$year) & depth > 20 & depth < 120) %>% nrow()
#> filter: removed 242,779 rows (97%), 6,674 rows remaining
pred_grid_sub <- filter(pred_grid, depth > 20 & depth < 120)
#> filter: removed 69,255 rows (28%), 180,198 rows remaining
pred_avg_sim_sub <- predict(mfull, newdata = pred_grid_sub, nsim = 100)
index_avg_sim_sub <- get_index_sims(pred_avg_sim_sub, area = rep(1/ncells_sub, nrow(pred_avg_sim_sub)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
# Combine and plot & compare
index_avg_sim_comp <- bind_rows(mutate(index_avg_sim, area = "full"),
mutate(index_avg_sim_sub, area = "subset"))
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> mutate: new variable 'area' (character) with one unique value and 0% NA
ggplot(index_avg_sim_comp, aes(year, log(est), ymin = log(lwr), ymax = log(upr), color = area)) +
geom_line() +
geom_ribbon(alpha = 0.4)
ggplot(index_avg_sim_comp, aes(year, est, ymin = lwr, ymax = upr, color = area)) +
geom_line() +
geom_ribbon(alpha = 0.4)
colnames(pred_grid)
#> [1] "X"
#> [2] "Y"
#> [3] "depth"
#> [4] "year"
#> [5] "oxy"
#> [6] "temp"
#> [7] "lon"
#> [8] "lat"
#> [9] "sub_div"
#> [10] "ices_rect"
#> [11] "density_saduria"
#> [12] "biomass_spr"
#> [13] "biomass_her"
#> [14] "biomass_spr_sd"
#> [15] "biomass_her_sd"
#> [16] "density_cod"
#> [17] "density_fle"
#> [18] "depth_rec"
#> [19] "temp_rec"
#> [20] "oxy_rec"
#> [21] "density_cod_rec"
#> [22] "density_fle_rec"
#> [23] "density_saduria_rec"
#> [24] "depth_sd"
#> [25] "temp_sd"
#> [26] "oxy_sd"
#> [27] "density_cod_sd"
#> [28] "density_fle_sd"
#> [29] "density_saduria_sd"
#> [30] "year_f"
#> [31] "median_biomass_sprat_across_sub_div"
#> [32] "median_biomass_herring_across_sub_div"
#> [33] "median_biomass_sprat_across_rect_in_sub_div"
#> [34] "median_biomass_herring_across_rect_in_sub_div"
#> [35] "n_rect"
#> [36] "ln_length_cm"
#> [37] "year_ct"
#> [38] "biomass_her_sc"
#> [39] "biomass_her_sd_sc"
#> [40] "biomass_spr_sc"
#> [41] "biomass_spr_sd_sc"
#> [42] "density_cod_sc"
#> [43] "density_cod_rec_sc"
#> [44] "density_fle_sc"
#> [45] "density_fle_rec_sc"
#> [46] "density_saduria_sc"
#> [47] "density_saduria_rec_sc"
#> [48] "depth_sc"
#> [49] "depth_rec_sc"
#> [50] "oxy_sc"
#> [51] "oxy_rec_sc"
#> [52] "temp_sc"
#> [53] "temp_rec_sc"
# Plot the sims
ggplot(index_avg_sims, aes(y = as.factor(year), x = .value)) +
geom_density_ridges_gradient(rel_min_height = 0.01, alpha = 0.8, color = "grey10", size = 0.2,
quantile_lines = TRUE, quantiles = 0.5) +
coord_flip()
#> Picking joint bandwidth of 0.00302
# index_avg_sims %>%
# filter(year %in% c(1993, 2019)) %>%
# ggplot(., aes(as.factor(year), .value)) +
# geom_violin(width = 1.3)
# Now reshape data again to calculate row-wise differences
year_diff <- index_avg_sims %>%
filter(year %in% c(1993, 2019)) %>%
#mutate(.value = log(.value)) %>%
pivot_wider(names_from = year, values_from = .value) %>%
mutate("2019-1993" = `2019` - `1993`) %>%
pivot_longer(2:4) %>%
ungroup()
#> filter: removed 2,500 rows (93%), 200 rows remaining
#> pivot_wider: reorganized (year, .value) into (1993, 2019) [was 200x3, now 100x3]
#> mutate: new variable '2019-1993' (double) with 100 unique values and 0% NA
#> pivot_longer: reorganized (1993, 2019, 2019-1993) into (name, value) [was 100x4, now 300x3]
#> ungroup: no grouping variables
# Plot these
p1 <- ggplot(year_diff) +
geom_density(data = filter(year_diff, !name == "2019-1993"),
aes(value, fill = name), alpha = 0.6) +
scale_fill_brewer(palette = "Dark2") +
coord_cartesian(expand = 0)
#> filter: removed 100 rows (33%), 200 rows remaining
p2 <- ggplot(year_diff) +
geom_density(data = filter(year_diff, name == "2019-1993"),
aes(value), alpha = 0.6, fill = "grey") +
geom_vline(xintercept = 0, linetype = 2)
#> filter: removed 200 rows (67%), 100 rows remaining
p1/p2
# Calculate quantiles for each level (to go in the results section)
year_diff %>%
group_by(name) %>%
summarise(median = quantile(value, probs = 0.5),
upr97.5 = quantile(value, probs = 0.975),
lwr2.5 = quantile(value, probs = 0.025))
#> group_by: one grouping variable (name)
#> summarise: now 3 rows and 4 columns, ungrouped
# Percent change in condition factor
year_diff %>%
group_by(name) %>%
pivot_wider(names_from = name, values_from = value) %>%
mutate(perc_change = ((`2019`-`1993`)/`1993`)*100) %>%
ungroup() %>%
summarise(lwr2.5 = quantile(perc_change, probs = 0.025),
median = quantile(perc_change, probs = 0.5),
upr97.5 = quantile(perc_change, probs = 0.975))
#> group_by: one grouping variable (name)
#> pivot_wider: reorganized (name, value) into (1993, 2019, 2019-1993) [was 300x3, now 100x4]
#> mutate: new variable 'perc_change' (double) with 100 unique values and 0% NA
#> ungroup: no grouping variables
#> summarise: now one row and 3 columns, ungrouped
# From these models, predict annual condition factor
# I will use the prediction grid with ALL covariates set to 0, incl. depth, temp and oxygen
# Grabbing the number of cells to help with calculating the average
ncells24 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 24) %>% nrow()
ncells25 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 25) %>% nrow()
ncells26 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 26) %>% nrow()
ncells27 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 27) %>% nrow()
ncells28 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 28) %>% nrow()
# Use the `area` argument here to turn the total into an average by giving it one over the number of cells
pred_avg_24_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 24), sims = 100)
#> Warning: The `sims` argument of `predict.sdmTMB()` is deprecated as of sdmTMB 0.0.21.
#> Please use the `nsim` argument instead.
pred_avg_25_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 25), sims = 100)
pred_avg_26_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 26), sims = 100)
pred_avg_27_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 27), sims = 100)
pred_avg_28_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 28), sims = 100)
index_avg_24_sim <- get_index_sims(pred_avg_24_sim, area = rep(1/ncells24, nrow(pred_avg_24_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_25_sim <- get_index_sims(pred_avg_25_sim, area = rep(1/ncells25, nrow(pred_avg_25_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_26_sim <- get_index_sims(pred_avg_26_sim, area = rep(1/ncells26, nrow(pred_avg_26_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_27_sim <- get_index_sims(pred_avg_27_sim, area = rep(1/ncells27, nrow(pred_avg_27_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_28_sim <- get_index_sims(pred_avg_28_sim, area = rep(1/ncells28, nrow(pred_avg_28_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
# Combine and plot
index_avg_24_sim <- index_avg_24_sim %>% mutate(sub_div = 24)
index_avg_25_sim <- index_avg_25_sim %>% mutate(sub_div = 25)
index_avg_26_sim <- index_avg_26_sim %>% mutate(sub_div = 26)
index_avg_27_sim <- index_avg_27_sim %>% mutate(sub_div = 27)
index_avg_28_sim <- index_avg_28_sim %>% mutate(sub_div = 28)
div_index <- bind_rows(index_avg_24_sim, index_avg_25_sim, index_avg_26_sim, index_avg_27_sim, index_avg_28_sim)
div_index2 <- bind_rows(mutate(div_index, sub_div = as.factor(sub_div)),
mutate(index_avg_sim, sub_div = as.factor("Total")))
pal <- c(brewer.pal(n = 5, name = "Dark2"), "grey10")
# Remember to log now that we use sims
cond_index <- ggplot(div_index2, aes(year, est, color = sub_div, linetype = sub_div)) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
labs(y = "Predicted Le Cren's condition factor", x = "Year") +
geom_line(alpha = 0.8) +
geom_point(alpha = 0) +
scale_color_manual(values = pal, name = "ICES subdivision") +
scale_linetype_manual(values = c(1,1,1,1,1,0)) +
geom_point(data = index_avg_sim, aes(year, est), inherit.aes = FALSE,
size = 2.5, shape = 21, fill = "gray1", color = "white") +
geom_errorbar(data = index_avg_sim, inherit.aes = FALSE,
aes(x = year, ymax = upr, ymin = lwr),
width = 0, alpha = 0.8, color = "gray1", size = 0.7) +
theme(axis.title = element_text(size = 9),
legend.position = "bottom",
legend.background = element_blank(),
legend.title = element_text(size = 9),
legend.spacing.x = unit(0, 'cm'),
plot.margin = unit(c(0.4, 0.4, 0.4, 0), "cm"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank()) +
guides(linetype = "none",
color = guide_legend(title.position = "top", title.hjust = 0.5, nrow = 1,
keywidth = 0.1, keyheight = 0.1, default.unit = "inch",
override.aes = list(linetype = c(1,1,1,1,1,0),
shape = c(NA,NA,NA,NA,NA,16),
alpha = rep(0.8, 6),
color = pal))) +
NULL
cond_index
# Ridge-plot of data
data_ridge <- ggplot(d, aes(x = le_cren, y = fct_rev(as.factor(year)), fill = factor(stat(quantile)))) +
stat_density_ridges(
rel_min_height = .01,
scale = 3,
geom = "density_ridges_gradient",
calc_ecdf = TRUE,
alpha = 0.1,
quantiles = c(0.25, 0.75)) +
scale_fill_manual(
name = "Probability",
values = c("grey50", "gray70", "grey90"),
labels = c("[0,0.25]", "[0.25,0.75]", "[0.75,1]")) +
xlim(0.65, 1.55) +
labs(x = "Le Cren's condidition factor",
y = "Year") +
guides(fill = guide_legend(nrow = 1, title.position = "top", title.hjust = 0.5)) +
theme(axis.title = element_text(size = 9),
legend.position = "bottom",
legend.key.size = unit(0.5, "line"),
legend.background = element_blank(),
legend.title = element_text(size = 9),
legend.spacing.x = unit(0.1, 'cm'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.margin = unit(c(0.4, 0.8, 0.4, 0), "cm")) +
NULL
data_ridge + cond_index + plot_annotation(tag_levels = 'A') +
plot_layout(widths = c(1, 2.2))
#> Warning: Removed 926 rows containing non-finite values (stat_density_ridges).
ggsave("figures/cond_index.pdf", width = 17, height = 17, units = "cm")
#> Warning: Removed 926 rows containing non-finite values (stat_density_ridges).
## Calculate delta_fulton and relate to delta_oxygen, by subdivision
div_index_delta_cond <- div_index2 %>%
group_by(sub_div) %>%
filter(year %in% c(1993, 1994, 1995, 2017, 2018, 2019)) %>%
dplyr::select(year, log_est, sub_div) %>%
pivot_wider(names_from = year, values_from = log_est) %>%
mutate(start_cond = mean(`1993`, `1994`, `1995`),
end_cond = mean(`2017`, `2018`, `2019`),
delta_cond = end_cond - start_cond) %>%
ungroup() %>%
dplyr::select(sub_div, delta_cond)
write.csv(div_index_delta_cond, "output/div_index_delta_cond.csv")
# What is one standard deviation of oxygen?
sd(d$oxy)
#> [1] 1.846734
# Extract random and fixed coefficients from the full model
mfull_est <- bind_rows(tidy(mfull, effects = "ran_par", conf.int = TRUE) %>%
filter(term %in% c("sigma_O", "sigma_E")),
tidy(mfull, effects = "fixed", conf.int = TRUE) %>%
filter(!grepl('year', term))) %>%
mutate(term = factor(term))
mfull_est <- mfull_est %>%
mutate(group = "Herring",
group = ifelse(term %in% c("biomass_spr_sc", "biomass_spr_sd_sc"), "Sprat", group),
group = ifelse(term %in% c("density_cod_rec_sc", "density_cod_sc"), "Cod", group),
group = ifelse(term %in% c("density_fle_rec_sc", "density_fle_sc"), "Flounder", group),
group = ifelse(term %in% c("density_saduria_rec_sc", "density_saduria_sc"), "Saduria", group),
group = ifelse(term %in% c("depth_rec_sc", "depth_sc"), "Depth", group),
group = ifelse(term %in% c("oxy_rec_sc", "oxy_sc"), "Oxygen", group),
group = ifelse(term %in% c("temp_rec_sc", "temp_sc"), "Temp", group),
group = ifelse(term == "sigma_E", "Random", group),
group = ifelse(term == "sigma_O", "Random", group)
)
mfull_est <- mfull_est %>%
mutate(scale = "Sub-division",
scale = ifelse(term %in% c("biomass_her_sc", "biomass_spr_sc", "density_cod_rec_sc", "density_fle_rec_sc",
"density_saduria_rec_sc", "depth_rec_sc", "oxy_rec_sc", "temp_rec_sc"), "Rectangle", scale),
scale = ifelse(term %in% c("density_cod_sc", "density_fle_sc", "density_saduria_sc",
"depth_sc", "oxy_sc", "temp_sc"), "Haul", scale),
scale = ifelse(term == "sigma_E", "Spatial/spatiotemporal s.d.", scale),
scale = ifelse(term == "sigma_O", "Spatial/spatiotemporal s.d.", scale)
)
# Plot effects
# This is the order changed labels should go in!
levels(mfull_est$term)
#> [1] "biomass_her_sc" "biomass_her_sd_sc" "biomass_spr_sc"
#> [4] "biomass_spr_sd_sc" "density_cod_rec_sc" "density_cod_sc"
#> [7] "density_fle_rec_sc" "density_fle_sc" "density_saduria_rec_sc"
#> [10] "density_saduria_sc" "depth_rec_sc" "depth_sc"
#> [13] "oxy_rec_sc" "oxy_sc" "sigma_E"
#> [16] "sigma_O" "temp_rec_sc" "temp_sc"
# I want the random effects to be dark gray...
pal <- brewer.pal(n = 8, name = "Dark2")
pal2 <- c(pal[1:5], "black", pal[6:8])
# Sort the terms so that random effects are at the top...
mfull_est <- mfull_est %>%
mutate(term2 = ifelse(term %in% c("sigma_E", "sigma_O"), 2, 1))
ggplot(mfull_est, aes(reorder(term, term2), estimate, color = group, fill = group, shape = scale)) +
geom_hline(yintercept = 0, linetype = 2, color = "gray40", alpha = 0.5) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
geom_point(size = 2.5) +
scale_color_manual(values = pal2, name = "") +
scale_shape_manual(values = c(15, 19, 21, 17), name = "") +
scale_fill_manual(values = rep("white", 9), name = "") +
labs(x = "", y = "Standardized coefficient") +
scale_x_discrete(breaks = levels(mfull_est$term),
labels = c(expression(Herring[rec]),
expression(Herring[sub]),
expression(Sprat[rec]),
expression(Sprat[sub]),
expression(Cod[rec]),
expression(Cod[haul]),
expression(Flounder[rec]),
expression(Flounder[haul]),
expression(Saduria~entomon[rec]),
expression(Saduria~entomon[haul]),
expression(Depth[rec]),
expression(Depth[haul]),
expression(Oxygen[rec]),
expression(Oxygen[haul]),
expression(sigma[E]),
expression(sigma[O]),
expression(Temp[rec]),
expression(Temp[haul])
)) +
coord_flip() +
theme(plot.margin = unit(c(0.4, 0.4, 0.4, 0), "cm"),
legend.position = "bottom",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank()) +
guides(color = "none", fill = "none", shape = guide_legend(ncol = 2)) +
NULL
ggsave("figures/effect_size.pdf", width = 17, height = 17, units = "cm")
# Just a test to see the labels were alright
#ggplot(mfull_est, aes(reorder(term, estimate), estimate)) +
# ggplot(mfull_est, aes(term, estimate)) +
# geom_hline(yintercept = 0, linetype = 2, color = "red", alpha = 0.5) +
# geom_point(size = 2) +
# labs(x = "", y = "Standardized coefficient") +
# coord_flip()
# Now save this for comparison in sensitivity analysis script
write.csv(mfull_est, "output/mfull_est.csv")
pred_map <- predict(mfull, newdata = pred_grid)
pal <- wes_palette("Zissou1", 21, type = "continuous")
plot_map_raster +
geom_raster(data = filter(pred_map, depth < 120 & depth > 10), aes(x = X*1000, y = Y*1000, fill = exp(est))) +
scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") +
facet_wrap(~ year, ncol = 5) +
theme_facet_map() +
NULL
ggsave("figures/supp/condition/all_years_condition_map_covars.png", width = 6.5, height = 6.5, dpi = 600)
# And a smaller map for selected years
plot_map_raster_labels +
geom_raster(data = filter(pred_map, year %in% c(1994, 2001, 2008, 2018) & depth < 120 & depth > 10),
aes(x = X*1000, y = Y*1000, fill = exp(est))) +
scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") +
facet_wrap(~ year, ncol = 2) +
theme_plot() +
NULL
ggsave("figures/condition_map_covars.pdf", width = 17, height = 17, units = "cm")
# Plot correlation between random effects and covariates
# Drop NA covariates for the correlation plot
# colnames(pred_map)
#
# pred_map_corrplot <- pred_map %>%
# drop_na(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
# density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc,
# density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
# oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc, omega_s, epsilon_st) %>%
# dplyr::select(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
# density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc,
# density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
# oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc, omega_s, epsilon_st)
#
# corr_re <- round(cor(pred_map_corrplot), 2)
#
# ggcorrplot(corr_re, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
# theme_classic() +
# theme(axis.text.x = element_text(angle = 90, vjust = 0.3))
#
# ggsave("figures/supp/condition/correlation_variables_random_effects.png", width = 6.5, height = 6.5, dpi = 600)
# Plot spatiotemporal random effect
plot_map_raster +
geom_raster(data = pred_map, aes(x = X * 1000, y = Y * 1000, fill = epsilon_st)) +
scale_fill_gradient2() +
facet_wrap(~ year, ncol = 5) +
ggtitle("Spatiotemporal random effects") +
NULL
ggsave("figures/supp/condition/epsilon_st_map.png", width = 6.5, height = 6.5, dpi = 600)
# Plot spatial random effect
plot_map_raster +
geom_raster(data = filter(pred_map, year == 1999), aes(x = X * 1000, y = Y * 1000, fill = omega_s)) +
scale_fill_gradient2() +
facet_wrap(~ year, ncol = 5) +
theme_plot() +
ggtitle("Spatial random effects") +
NULL
ggsave("figures/supp/condition/omega_s_map.png", width = 6.5, height = 6.5, dpi = 600)
# Fit a linear model to each prediction grid of the estimate over time
# https://community.rstudio.com/t/extract-slopes-by-group-broom-dplyr/2751/7
time_slopes_by_year <- pred_map %>%
mutate(id = paste(X, Y, sep = "_")) %>%
split(.$id) %>%
purrr::map(~lm(est ~ year, data = .x)) %>%
purrr::map_df(broom::tidy, .id = 'id') %>%
filter(term == 'year')
# Plot the slopes
time_slopes_by_year2 <- time_slopes_by_year %>%
separate(id, c("X", "Y"), sep = "_") %>%
mutate(X = as.numeric(X),
Y = as.numeric(Y))
plot_map_raster +
geom_raster(data = time_slopes_by_year2, aes(x = X * 1000, y = Y * 1000, fill = estimate)) +
scale_fill_gradient2(midpoint = 0) +
facet_wrap(~ term, ncol = 5) +
theme_plot() +
ggtitle("Time slopes by each pred grid") +
NULL
<img src=“condition_model_cf_files/figure-html/calculate”spatial trend”-1.png” width=“672” style=“display: block; margin: auto;” />
ggsave("figures/supp/condition/spatial_trend_condition.png", width = 6.5, height = 6.5, dpi = 600)
# Prepare prediction data frame
nd_dep <- data.frame(depth_sc = seq(min(d$depth_sc), max(d$depth_sc), length.out = 100))
nd_dep <- nd_dep %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model (AIC-selected)
p_margin_dep <- predict(mfull, newdata = nd_dep, se_fit = TRUE, re_form = NA)
mar_dep <- ggplot(p_margin_dep, aes(depth_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
theme(aspect.ratio = 1) +
labs(x = expression(Depth[haul]))
# Prepare prediction data frame
nd_oxy <- data.frame(oxy_sc = seq(min(d$oxy_sc), max(d$oxy_sc), length.out = 100))
nd_oxy <- nd_oxy %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
#oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model (AIC-selected)
p_margin_oxy <- predict(mfull, newdata = nd_oxy, se_fit = TRUE, re_form = NA)
mar_oxy <- ggplot(p_margin_oxy, aes(oxy_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se)))+
geom_ribbon(alpha = 0.4) +
geom_line() +
theme(aspect.ratio = 1) +
labs(x = expression(Oxygen[haul]))
# Prepare prediction data frame
nd_tem <- data.frame(temp_sc = seq(min(d$temp_sc), max(d$temp_sc), length.out = 100))
nd_tem <- nd_tem %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
#temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model (AIC-selected)
p_margin_tem <- predict(mfull, newdata = nd_tem, se_fit = TRUE, re_form = NA)
mar_temp <- ggplot(p_margin_tem, aes(temp_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
theme(aspect.ratio = 1) +
labs(x = expression(Temperature[haul]))
# Prepare prediction data frame
nd_cod <- data.frame(density_cod_sc = seq(min(d$density_cod_sc), max(d$density_cod_sc), length.out = 100))
nd_cod <- nd_cod %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
#density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model (AIC-selected)
p_margin_cod <- predict(mfull, newdata = nd_cod, se_fit = TRUE, re_form = NA)
mar_cod <- ggplot(p_margin_cod, aes(density_cod_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
theme(aspect.ratio = 1) +
labs(x = expression(Cod[haul]))
# Prepare prediction data frame
nd_spr <- data.frame(biomass_spr_sd_sc = seq(min(d$biomass_spr_sd_sc), max(d$biomass_spr_sd_sc), length.out = 100))
nd_spr <- nd_spr %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
#biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model (AIC-selected)
p_margin_spr <- predict(mfull, newdata = nd_spr, se_fit = TRUE, re_form = NA)
mar_spr <- ggplot(p_margin_spr, aes(biomass_spr_sd_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
#coord_cartesian(ylim = c(-4.65, -4.4)) +
theme(aspect.ratio = 1) +
labs(x = expression(Sprat[sub]))
# Prepare prediction data frame
nd_sad <- data.frame(density_saduria_sc = seq(min(d$density_saduria_sc), max(d$density_saduria_sc), length.out = 100))
nd_sad <- nd_sad %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
#density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model (AIC-selected)
p_margin_sad <- predict(mfull, newdata = nd_sad, se_fit = TRUE, re_form = NA)
mar_sad <- ggplot(p_margin_sad, aes(density_saduria_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
#coord_cartesian(ylim = c(-4.75, -4.4)) +
theme(aspect.ratio = 1) +
labs(x = expression(Saduria~entomon[haul]))
# mar_dep + mar_oxy + mar_temp + mar_cod + mar_spr + mar_sad
# plot_annotation(tag_levels = 'A')
#
# ggsave("figures/supp/condition/marginal_effects.png", width = 6.5, height = 6.5, dpi = 600)
p_margin_dep2 <- p_margin_dep %>%
mutate(variable = "Depth (haul)") %>%
dplyr::select(est, est_se, depth_sc, variable) %>%
rename(value = depth_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_margin_oxy2 <- p_margin_oxy %>%
mutate(variable = "Oxygen (haul)") %>%
dplyr::select(est, est_se, oxy_sc, variable) %>%
rename(value = oxy_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_margin_tem2 <- p_margin_tem %>%
mutate(variable = "Temperature (haul)") %>%
dplyr::select(est, est_se, temp_sc, variable) %>%
rename(value = temp_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_margin_cod2 <- p_margin_cod %>%
mutate(variable = "Cod (haul)") %>%
dplyr::select(est, est_se, density_cod_sc, variable) %>%
rename(value = density_cod_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_margin_spr2 <- p_margin_spr %>%
mutate(variable = "Sprat (sub)") %>%
dplyr::select(est, est_se, biomass_spr_sd_sc, variable) %>%
rename(value = biomass_spr_sd_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_margin_sad2 <- p_margin_sad %>%
mutate(variable = "S. entomon (haul)") %>%
dplyr::select(est, est_se, density_saduria_sc, variable) %>%
rename(value = density_saduria_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
margs <- bind_rows(p_margin_dep2, p_margin_oxy2, p_margin_tem2,
p_margin_cod2, p_margin_spr2, p_margin_sad2)
ggplot(margs, aes(value, exp(est), ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
facet_wrap(~variable, scales = "free_x") +
labs(y = "Le Cren's condidition factor",
x = "Scaled value") +
theme_plot() +
theme(axis.text.x = element_text(angle = 0)) +
NULL
ggsave("figures/supp/condition/marginal_effects.png", width = 6.5, height = 6.5, dpi = 600)